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Thermodynamic engines are focusing increasing attention in the context of solar-based electric power 
generation. Knowledge-based models of such engines are sometimes difficult to derive and when they 
are available, their simulation may be numerically a rather heavy task given the control updating period 
that may be needed. In the present work a generic nonlinear identification framework that enables the 
dynamics of the key quantities of a thermodynamic engine to be captured is proposed. Such a fast model 
can then be used in the simulation and the control design stage of the whole electric power generation 
station. The proposed identification framework is validated on a recently developed knowledge-based 
model of a beta-type Stirling engine with rhomic-drive mechanism. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Thermodynamic engines transform the thermal energy into 
mechanical energy. This is done through a positive energy cycle 
in which a fluid circulates between a hot source and a cold source. 
The resulting mechanical energy can therefore be used to drive an 
alternator which is followed by a voltage conditioning stage leading 
to a usable electrical energy by the end user loads. 

A typical view of a thermodynamic engine in such a global hybrid 
power station framework is shown in Fig. 1. 

Running such hybrid power stations is today an important chal¬ 
lenge as it can be attested by the recent works involving the use of 
renewable wind/photovoltaic related energy sources [7,1,10,11]. 
One reason is that more than 1.4 billion people have no access 
to electricity, mainly in Asia and Africa. Classical alternative fuel- 
based solutions such as diesel engines suffer from many drawbacks 
such as their low reactivity to load changes, their high emissions 
and noise to mention only few issues [8]. 

Roughly speaking, a thermodynamic engine in a hybrid power 
generation station provides the medium term power adaptation 
while an electric storage element (super-capacitor and/or battery) 
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accommodates for sudden changes in the load power demand. It 
results (see [4 for a detailed demonstration) that the dynamics of 
the thermodynamic engine highly determine the size and the type 
of the storage element. It determines also to which extent one can 
avoid the use of batteries which is becoming nowadays necessary 
for obvious ecological reasons. 

Now in order to perform such a closed-loop quantitative analysis 
on a variety of scenario, it is necessary to dispose of a fast simula¬ 
tor of the thermodynamic engine in its loop. This can be obtained 
either using a knowledge-based model or using a black-box model. 
Note also that even when a knowledge-based model is available, 
it can be used to derive black-box model in order to accelerate the 
simulation and the control algorithm development task. 

To the authors knowledge, several companies are currently 
developing thermodynamic engines following the basic concepts of 
positive energy cycles (Stirling, Rankine, etc.). Based on some ele¬ 
mentary computation schemes, the main parameters of the engine 
can be computed given the targeted output mechanical power 
range. It comes out that thermodynamic engines can be available 
before a dynamic model is derived. In this case, real-life experi¬ 
ments can be conducted to derive black-box models without the 
need for a mathematical knowledge-based model. Note however 
that knowledge-based models are mandatory for process optimiza¬ 
tion in dynamic context. Black-box models can only reproduce a 
predesigned process behavior. 

The aim of this paper is to give a generic easily identifiable non¬ 
linear structure that captures the dynamics of the key quantities of 
a thermodynamic engine which are necessary to simulate it in its 
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Fig. 1 . An example of a hybrid power station with its thermodynamic engine (here a beta-type Stirling engine) followed by the electric power conditioning stage. 


environment. The adopted structure has been recently proposed in 
[2,3] where it has been used to derive a fast moving-horizon state 
observer [2] and to derive a dynamic model of NOx emission in a 
diesel engine [3 . 

In the absence of experimental data, the proposed identifica¬ 
tion scheme is validated in the present contribution using the very 
detailed model of a beta-type Stirling engine recently proposed in 
[5]. This knowledge-based model takes into account, among other 
non-ideal features, the non-isothermal effects, the effectiveness of 
the regenerative channel, and the thermal resistance of the heating 
head. This is why such a detailed model can reasonably emu¬ 
late an experimental data generator for the identification method 
proposed in this paper. This holds as far as only the presumably 
measured signals are used to build the black-box model. 

It is worth underlying that although the above cited specific 
Stirling engine is used throughout the paper, the identification 
structure and methodology are quite general and should be suc¬ 
cessfully used to derive a databased model for any thermodynamic 
engine. This is because almost all thermodynamic engines involve 
the same components and describe very similar cycles in the (P, V) 
coordinates. For instance, the framework proposed in the present 
paper has been also successful in identifying a Rankine-cycle-based 
engine. However, the corresponding results cannot be published for 
industrial reasons. 

The contributions of the papers can be summarized by the fol¬ 
lowing items: 

(1) The paper suggests a systematic parsimonious identification 
problem setting that targets only those variables that are rele¬ 
vant for power stations simulation and control. In the specific 
case of the Stirling engine considered in the present paper, 
these variables are the mean work and the mean heat exchange 
during an engine cycle. Note that this problem definition is 
independent of the method and the identification structure that 
can be used to solve the resulting identification problem. 

(2) The second contribution of the paper is to show that using a 
recently developed identification framework, the identification 
problem described in the preceding item leads to astonishingly 
faithful models of the dynamics of the targeted key quantities 
(work and heat exchange) that need to be captured. It is also 
conjectured (for the reasons cited above) that this success is 
likely to be encountered on a whole family of thermodynamic 
engines that are currently widely developed for green energy 
generation and conversion. 

It is by no means the claim of the present paper that the identifi¬ 
cation framework used herein is the only one capable of producing 


nice matching results. It is true however that the proposed frame¬ 
work does show nice properties that are recalled in Section 3. 

The paper is organized as follows: 

• In Section 2, the identification problem that has to be solved in 
order to derive a parsimonious dynamic model of the thermody¬ 
namic engine is stated. By parsimonious, it is meant a dynamic 
model that incorporates only the necessary elements that are 
needed to interact with the remaining blocks of the hybrid power 
generation station (see Fig. 1) and to design the feedback laws. 

• In Section 3, the general identification scheme proposed in 2,3] 
is briefly recalled and its application in the specific identification 
context described in Section 2 is discussed. It is then shown that 
the identification problem can be put into a general nonlinear 
static map identification. 

• Section 4 gives a brief sketch of the knowledge-based model of 
the beta-type Stirling engine with rhomic-drive mechanism that 
has been recently proposed in 5] and used in the next section as 
a data generator for the identification task. 

• The identification framework is validated in Section 5 where 
the learning/validation data are described and the identification 
results are shown. 

• Finally, Section 6 concludes the paper and give hints for further 
investigations. 

2. The identification problem 

In this section, the necessary elements that need to be identified 
regarding the thermodynamic engine in order to build the simula¬ 
tor of the whole hybrid power generation station are analyzed. This 
enables the identification problem to be clearly stated. 

The hybrid power station model can be shortly written as fol¬ 
lows: 

IQ - r e {£2, x e , u e ) + r m {£2,x th , u th , w th ) - /}(...) (1) 

x e = F(x e , £2, u e , w e ) (2) 

where Q (rad/s) is the angular velocity of the main shaft that is 
shared by the engine and the alternator. x e and u e are the state 
and the control vectors of the electrical part. x th and u th are the 
state and the control vectors of the thermodynamic engine. r e is 
the torque applied on the shaft from the electric part while r m is 
the mechanical torque that results from the Stirling engine posi¬ 
tive energy cycle. Fj is the friction torque. Finally, the vectors w e 
and w th represent the non controlled signals acting on the electri¬ 
cal part (user power, disturbances) and the thermodynamic part 
(hot temperature T H for instance). For the complete definition of F, 
interested readers can refer to the description given in [4 . 
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A typical instantiation of the electrical equation (2) can be 
found in [4 where a voltage conditioning stage is used including a 
super-capacitor. In the present paper, this part of the power gener¬ 
ation station is not considered since the interest is focused on the 
derivation of a fast and parsimonious model of the thermodynamic 
engine. 

It comes out that the rotational speed £2 is the coupling quantity 
between the mechanical and the electrical parts and that the only 
signal that summarizes the behavior of the thermodynamic engine 
is the torque r m that can be linked to the power P m delivered by 
the engine through the following equation: 


r m (Q,x th , 



P m (Q, x th , u th , w th ) 
Q 



Now given the mechanical inertia, the mean power P m during a 
cycle can be used instead of the instantaneous power. But the mean 
power P m can be given by: 



. _ 2n 



where W ou t is the mechanical work produced during a single ther¬ 
modynamic cycle. This leads to: 



1 WoMx th ,u th ,w th ) 


2n 



where the arguments are the mean values of the original arguments 
over the current cycle. 

This means that as far as simulation of the whole system is con¬ 
cerned, the key thermodynamic quantity to be identified is the 
work produced by the thermodynamic engine during a cycle for 
a given mean rotational speed, given mean state and given mean 
input vectors. 

As mentioned before, one of the reasons for which a fast dynamic 
model is needed lies in the control design issue. Indeed, it is worth 
operating the engine in a high efficiency region in the (power, 
speed) map. This map generally shows a bell-like curve that may 
depend on the input u th for a given rotational speed £2. 

Therefore, beside the identification of the map W out involved in 
(5) one needs to identify the efficiency map: 


W 0U t(C2, x th , u th , w th ) 

r l = - — _ f h -th —th ( 6 ) 

Qi n (n, X th , U th , W th ) 

where Qi n is the heat exchange with the hot source during a cycle. 

The expressions of W ou t and Qj n in terms of the measured 
quantities depend on the specific engine at hand. In Section 4, the 
expressions of W ou t and Qin (and hence the efficiency rj) are given for 
the specific case of the beta-type Stirling engine with rhomic-drive 
mechanism. 


2.2. Notation 


In the remainder of this section and in order to simplify the 
expressions, the vectors u is temporarily used to denote both the 
rotational speed, the controlled and the exogenous signal of the 
thermodynamic engine, namely: 


/«\ 

u th 

y w th J 



The rationale behind this is that from the identification point of 
view, there is no need to distinguish these inputs. Note that during 
the identification experiment, the rotational speed can be imposed 
using controlled motor as a load and can therefore be viewed as 
exogenous signal. Moreover, the superscript th is omitted since only 
the thermodynamic engine’s variables are concerned in the sequel. 


Therefore, the two quantities to be identified are given by: 


W ou t(x, u); Qin(x, u ) 



since then r)(x, u) is obtained through (6). Now since the two 
quantities represent the produced mechanical work and the heat 
exchange evaluated on a cycle, cycle-related notation is necessary 
to properly express the identification problem. These are given 
below. 

Note first of all that a cycle duration is defined to be the time 
interval between two successive instants t k and t fc+1 such that: 


6{t k ) e {2m i e N}; 0(t k+ j) = 6{t k ) + 2n 



where 0 is the angle of the Stirling engine shaft. 

Let us denote by t^ the successive measurement acquisition 
instants inside the time interval [t k , t k +\ ], namely: 


k < tf < r[ !+i) < 


< t k k < t k+ , 


( 10 ) 


Note that the number N™ of such instants does depend on k since 
the length of a cycle is not constant and depends on the rotational 
speed profile £?(•) around the instant t k . More precisely: 


* f /<+1 


2tt = 


£2(t)dt 


( 11 ) 


Given a cycle number /<, the following notation: 


x(t, /<); u(t,k ); £2{t,k)- te[t k ,t k+1 ] (12) 

is used hereafter to refer to the value of the state, control and rota¬ 
tional speed over the cycle duration. 

The mechanical work and the heat exchange during the cycle 
number k are denoted by W ou t{k) and Qm(/<). One obviously has: 


W out (k) = W(x(t /o /<), u( •, /<)); Qj n (fc) = Q(x(t k , /<), u( •, /<)) (13) 


since everything that happens during a cycle depends on the initial 
state x(t k , k) at the beginning of the cycle and the input profile u(-, 
k) during the cycle. 


2.2. Removing the state vector 

One of the obviously difficult issues in identifying the maps (8) 
lies in the presence of the state vector x. Indeed, if the entire state 
vector is used in the identified map, then the use of the black-box 
model would need the whole state vector to be reconstructed which 
may be cumbersome. Fortunately, this difficulty can be overcome 
thanks to the inherent open-loop stability of the thermodynamic 
engines. 

Indeed, open-loop stability of cyclic processes implies that the 
effects of initial state disappear asymptotically as time progresses. 
Using the notation above, this property can be formally written as 
follows: 

x(t fc ,k) = G(u(.,k-N + l),...,u(.,fc-l)) + o(l) (14) 

which can be read informally as follows: 

The value of the state at the beginning of a cycle can be approxi¬ 
mately viewed as a function of the input profiles during the last 
N cycles provided that N is sufficiently high. 

Now using (13) together with (14) enables to write 


Wout(k) « W(C(u( •, k — N + 1),..u( •, k — 1)), u( •, k)) (15) 

Qin(k) « Q(G(u( •, k - N + 1),..u( •, k- 1)), u( •, /<)) (16) 

Now assuming that the variations of the components of the input 
u are slow when compared to the duration of the cycle, one can 
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replace u(■ , k) by the mean value u(k) during the period, this leads 
to: 


WoutW « W(G(u(k - N + 1),.. u(k - 1)), u(k)) (17) 

Qin(fe) ~ Q(G(u(k - N + 1),.. u(k - 1)), u(k)) (18) 

which can be written in a more compact form using obvious defi¬ 
nitions of Fi and F 2 : 


W out (/<) = F-i (Z(l<)); Qj n (/<) = F 2 (Z(/<)) 


where Z(/<) is defined by: 


Z(k): 


/ u(k) \ 

: el nz ; 

\u(k-N + l) / 


n z :=N n u 


(19) 


( 20 ) 


In the sequel, the vector Z is referred to as the regressor. Anticipat¬ 
ing on the presentation of Sections 4 and 5, the components of the 
vector u involved in the definition of the regressor Z are given by: 


u := 



( 21 ) 


where T H (K) is the temperature of the hot source and <5 is a 
parameter that controls the flow rate of the air in the regenera¬ 
tor. Obviously, one can add also the cold temperature T L but this is 
considered here to be constant as it does not depend on the sunlight. 

In the following section, a recently proposed general framework 
[3] for the identification of nonlinear maps such as F\ and F 2 is 
recalled. 


3. Identification of general nonlinear relationships 

The identification of the two maps F\ and F 2 involved in (19) 
gathers two instantiations of the following general nonlinear iden¬ 
tification problem: 

Definition 1 (The Identification Problem). Given the data: 

Mk), Z(k))} N k lT (22) 


Note that since r is strictly monotonic, does exist. The ratio¬ 
nale behind this restriction is that the original problem of finding F 
can be replaced by the one consisting in finding L e M nz and r such 
that: 

r{q)^Z T l (25) 

Now assume that the map r is parametrized using the following 
structure: 

n b 

nq) = y[EfiK^q))lnf. m := g ~ g ™ in (26) 

z —' 4 max — 4 min 

J=1 

in which q min and q max are the lower and the upper values of q over 
the identification data (22) while B^(-) are a set of functions that 
are defined on [0,1]. 

The detailed analytical definition of the maps B^(-) is given in 
Appendix A.l. Fig. 2 shows the allure of these maps for = 0.5. 

Remark 1. The choice of the function basis shown in Fig. 2 is not 
unique. Nevertheless, this choice is motivated by the monotonicity 
assumption on FT Indeed, any one of the curves shown in Fig. 2 
is a monotonic function that can be viewed alone as a candidate 
function for r while each one of them would need many terms 
of standard basis functions such as polynomials to be correctly 
represented. 

Note that the expression of r(q) can be put in a matrix form 
using the vector of parameter ft e R nb and an obvious definition of 
m e M lxn fc as follows: 

r(q) = m(q))].fi (27) 

which enables (25) to be transformed into: 

m(q)),-Z T ]- J^O (28) 

that has to be solved in least squares sense under the following two 
linear constraints on the vector of coefficients ft 

• The first constraint expresses the fact that r has to be strictly 
monotonic: 


find a map F : M nz -> M that approximately maps Z to q in a least 
squares sense, that is: 


Vf €[0,1], 



■ ft > € 


(29) 


q « F(Z) (23) 

Identification of nonlinear relationships is an open issue. 
Many frameworks have been proposed including neural networks, 
Wiener-Fiammerstein, Volterra series based formulations [9,6] to 
cite but few possibilities. 

Although any of these methods can be used to solve the identi¬ 
fication problem at hand, the method used in the sequel is the one 
that has been recently proposed in [3] since it shows the advantage 
of being based on the solution of a Quadratic Programming (QP) 
problem. The consequences of this feature are the following: 

(1) Since the solution of QP problems can be done extremely fast, 
it is possible to perform on-line updating of the model during 
the system life-time. 

(2) QP problem solution does not suffer from possible local min¬ 
ima. This feature can be used in an outer-loop in which the 
parameters of the identification structure can be optimized as 
it is explained later in this section. 

The price to pay in order to put the problem in a QP form is to 
restrict the class of functions F to the one given by: 

F(Z) = F" _1 (Z r L); r( •) strictly increasing 


• The second constraint is a normalization constraint that avoids 
the trivial solution ft = 0, L = 0: 

p l 

[J B(g)dg] • At = qmin j gmax (30) 





(24) 


Fig. 2. Allure of the maps fP( •) for /3 = 0.5 and n m = 10. 
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Note that (29) and (30) can be transformed into a finite number 
of linear constraints by defining a sufficiently dense grid of values 
of § over the interval [0, 1 ]. 

To summarize, the two optimal parameters of the solution, 
namely! e M nz and \X e R n b are obtained by solving the constrained 
linear least squares defined by (28)-(30). 

Note that the degree of nonlinearity of the solution can be eval¬ 
uated by the indicator: 

max ?s[o,i][(dB)/(d£X£)] - min ?€( 0 il] [(dB)/(d£)(£)] 
y « : - ma % | 0 ,||[ d »)/(^)(j)| - (3 » 

which is obviously a number that belongs to [0, 1 ] and that equals 
0 for linear models. Note that the denominator in (31) is obviously 
non vanishing thanks to the monotonicity constraint (29). 

Remark 2. Regarding the choice of the parameters N and n b 
involved respectively in the definition of the regressor size and the 
number of functions used to represent the map r, the following 
procedure can be used: 

(1) Choose maximum values n™ ax and N max on these parameters. 

(2) Solve the QP problem on the grid of possible values defined by 
the set: 



Fig. 3. Schematic view of the beta-type Stirling engine modeled in [5] and which is 
used in this paper to validate the proposed identification framework. 


{l,...,n^ ax }x{l,...,N max } 

This can be done extremely rapidly since the optimization prob¬ 
lem is a simple QP. 

(3) Denote by r min the value of the lowest residual which is typically 
obtained for the pair (n™ ax x N max ). 

(4) Fix a threshold 6 > 0. 

(5) Compute the pair (N, n m ) that is minimal (in lexicographic 
sense) which corresponds to a residual that is not greater than 
r min by more than 6%. The use of lexicographic order enhances 
smaller values of N since high values of N increase the regressor 
dimension. 


4. Description of the knowledge-based data generator 

As it is explained above, the formulation of the identification 
problem and the proposed identification framework are validated 
using data issued from a detailed knowledge-based model of a beta- 
type Stirling engine. For a detailed discussion of the model, the 
reader is invited to examine the recent paper 5] where the com¬ 
plete derivation of the model is given. Here, only a brief description 
of the model derivation is given in order to underline some specific 
features. 

This Stirling engine contains two compartments which are the 
expansion chamber and the compression chamber (see Fig. 3). 
These chambers are respectively close to the hot source with tem¬ 
perature T h and the cold source with a temperature T L . The air 
circulates between these two chambers through the regenerative 
channel. The volume of the chambers is uniquely determined by 
the instantaneous value of the angle 0 and the fact that the lengths 
i c and i e are constant. 

The model proposed in [5] assumes that the air obeys the ideal 
gas equation, which means that only two of the three variables (T, 
P, m) in each chamber are independent for a given value of 0 (and 
hence the volume). 

The authors of [5] choose the pair (T, m) in each chamber as the 
pair of independent variable leading to a state x th of dimension 5 
that is defined by: 


(32) 


while the exogenous input vector u th is given by: 
u := {£2 T h <5) r (33) 


Remark 3. In operational Stirling engines, the flow rate through 
the regenerative channel between the expansion and the compres¬ 
sion chambers is monitored through valves in order to control the 
heat transfer. In the simplified scheme of the Stirling engine con¬ 
sidered in the present paper the variable <5 seems to be constant by 
construction (see Fig. 3). However, once the model is obtained, one 
can emulate the effects of controlled flow rate on the efficiency of 
the engine by simulating dynamic variations of <5 in order to check 
whether the identification structure would be able to capture the 
resulting dynamics. 


The equations of the model take the following form: 


dll 

. 


dm ( 

dt ^ 

= Qin,e ~ 

H/out,e - 

dt X ( 

dU 




dt lc 

= Qin,c ~ 

W 0 ut,c ~ 

dt X ( 

dm 

dm 

dm 

dm 

dt ^ 

dt ’ 

dt ' L 

~ + dr 


v 

~2 


* 4 


d0 


= £2 


(34) 

(35) 

(36) 


where U e and U c represent the internal energies of the air in each 
chamber, Q n?e = {J H - T e )/R tl and Q n?c = (T L - T c )/Rt 2 are the heat 
exchanged by the chamber with the hot and the cold sources 
respectively. R tl and R t2 are the thermal resistances of the heating 
head and the cooling jacket respectively. W e?0 ut and W c?ou t are the 
mechanical works developed in the expansion and the compression 
chambers respectively. Their expressions involve the pressures P e 
and P c and the variations of the volumes V e and V c that depend 
exclusively on the current position 0 and its derivative £2. 

The dm\dt is the mass flow rate through the regenerative chan¬ 
nel (see Fig. 3). The pairs (h, v)\ e and (h, v)\ c represent the enthalpy 
and the velocity of the air at the boundaries between the chambers 
and the regenerative channel. The velocities v e and v c are linked to 
the mass flow rate dm\dt through the expressions: 


|(dm)/(dt)|. = |(dm)/(dt)| 

2p e 7ir8 ’ c 2p c 7Tr8 


(37) 


with r being the radius of the cylinder (see Fig. 3). 


x th := (Je T c m e m c 0) T 
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The expression of the mass flow rate (dm)/(dt) is obtained by 
writing the equation of momentum of a fully developed annular 
flow of an incompressible fluid. This means that the velocity is 
zero at the boundaries while it depends on the radial coordinate 
in between (see Fig. 3). The expression involved the instantaneous 
pressure difference between the two chambers and the velocity in 
the displacer (see 5] for the detailed expressions). 

Now expressing the internal energies in terms of the tem¬ 
peratures U = C V (T)-T, injecting it into Eqs. (34) and (35) and 
performing implicit discrete integration leads to equations of the 
form (in which x + and p + denote the values of x and p at the next 
sampling instant). 

x + =f(x, u,p+); p:=(P e ,P c ) (38) 

P = g6) (39) 

where (39) expresses the ideal gas equation in the expansion and 
the compression chambers. This obviously defines a fixed-point 
iteration in the pressures p = (P e , Pc) that has to converge for each 
time step of the model integration: 

P (i+1) =gof(x, u,p (0 ) =: C M (p (i) ) (40) 

As it has been mentioned in [5 , the convergence of the simulator 
needs extremely small sampling periods, our experience suggests 
that the sampling period: 

t = 5 x 10 -6 s 

seems to be an upper bound. This is what makes the simulator quite 
heavy for potential integration in a whole simulation framework of 
the hybrid power generation station. 

Fig. 4 shows a typical cycle of the beta-type Stirling engine 
machine described above. The cycle clearly shows how far it is from 
the ideal Stirling cycle which is composed of ideal isotherms and 
isobar curves. The cycles clearly show the positive sign of the out¬ 
put work as the clockwise cycle shows greater surface than the 
anti-clockwise cycle. 

( t)dt + I P c (t)V c (t)dt (41) 

Jo 



Typical cycle in the P-V plan 



Fig. 4. Typical (P-V)-c ycle of the beta-type Stirling engine used as a data generator 
for the validation of the proposed identification scheme. 


As for the heat exchanged with the hot source which is involved 
in the definition of the efficiency (6), it is given by the following 
expression as far as this Stirling engine is concerned: 




Rt\ 


(42) 


which means that given the geometry of the machine that enables 
the volume to be tracked as a function of 0 , the computation of W ou t 
and Qin needs the measurement of the pressures P e and P c as well 
as the temperature T e . 


5. Validation of the identification framework 

The identification framework recalled in Section 3 is used here 
to identify two maps F\ and F 2 such that: 

Qi ~ ^t(Z); q 2 =F 2 (Z) (43) 

where q\ = W ou t and q 2 = Qin while the regressor Z is given at each 
instant k by (20) in which u stands for the mean value of the vector 
u given by (21). In the following sections, the generation of the 
learning data (22) and the choice of the identification parameters 
are shown, the identification results are given before the identified 
solution is validated using validation data. Finally, the computation 
times of the knowledge-based model and the black-box model are 
compared. 

5.1. The learning data 

The learning data is generated by simulating the knowledge- 
based model described in Section 4 during 40 s using the input 
profiles given in Fig. 5. This leads to a scenario that contains 
N data = 646 cycles that are then used to build the regressor instances 
Z(/<) as well as the corresponding instances of q-i(k) = W 0Ut (k) and 
Q2W = Qin{k). 

Note that the learning data is based on the following bounds on 
the input signals: 

F2 e [40,120] (rad/s); 8 e [0.15, 0.25](mm); T H e [600, 900](K) 


5.2. The identification parameters 

The following parameters have been used in the nonlinear 
parameter structure: 

n m = 3; N = 3; 6 = 10 (44) 

which leads to a regressor length 
n z = N ■ n u = 9 

and a total number of model parameters equal to 
n p = n b + n z = 2 n m + n z = 15 parameters 

when reported to the number of cycles, this gives a compression 
factor of: 

472 

Compression factor = -j-^->43 (45) 

Note that the value N= 3 means that the prediction of the work 
Wout and the heat exchange Qj n during a cycle depends on the mean 
values of the input during the last three cycles including the current 
one. Recall that this feature is intimately linked to the open-loop 
stability of the engine as it has been explained in Section 2. 
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Fig. 5. The time profiles of the exogenous signals that are used to build the learning scenario. 
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Fig. 6. Result of the identification on the learning data. Note that while the map W out = F i (Z) is almost affine, higher nonlinearity is needed to approximate heat exchange 


Qln=F 2 (Z). 


5.3. Identification results on the learning data 

The identification results are shown in Fig. 6. More precisely the 
top subplots show the results regarding the identification of the 
work while the bottom subplots shows the results for the identi¬ 
fication of the efficiency. On each row, three subplots are given: 
at the left, the comparison between the true and the identified 
values. At the center, the allure of the map T(§) is given while 


the righter plot shows the derivative w.r.t. § of the same func¬ 
tion in order to easily appreciate the degree of nonlinearity of the 
function. 

Fig. 7 shows a temporal comparison of the true and the identi¬ 
fied quantities. Fig. 8 shows the comparison between the true and 
the identified efficiency. The latter is computed according to (6) in 
which the identified work W ou t and the identified heat exchange 
Qin are used. 


True and identified work 


True and identified Heat Exchange 



Fig. 7. Identification on the learning data: time profiles of the true and the identified work and heat exchange. 
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True and identified efficiency 



Fig. 8. Identification on the learning data: time profiles of the true and the identified 
efficiency. 

5.4. Validation of the black-box model 

As in any identification task, the validation of the black-box 
model must be done on temporal profiles of the exogenous sig¬ 
nals that differ from the ones that are used in the learning data. 
Fig. 9 shows the exogenous input profiles that define the validation 
scenario. 

The predicted profiles of the work and the heat exchange are 
compared to the true ones in Fig. 10 while the comparison between 
the predicted efficiency and the true one is shown in Fig. 11. 

The results clearly show the quality of the black-box model as 
it captures the dynamics of the identified variables in a totally dif¬ 
ferent scenario. 


true and identified heat exchanged 



Fig. 11. Identification on the learning data: time profiles of the true and the pre¬ 
dicted efficiency. 

Table 1 

Comparison of the computation time needed to simulate the scenario of Fig. 5 with 
the knowledge-based model of [5] and the black-box model. The knowledge-based 
model has been compiled in C language. Platform: Mac OSX, 2.3 GHz intel Core i7, 
8Go 1333 MHz DDR3. 

Computation times (scenario of Fig. 5) 

Black-box model 0.07 s 

Knowledge-based model 39 s 


has been compiled in C language in order to accelerate the compu¬ 
tation time. The comparison suggests an acceleration factor of 570. 

6. Conclusion and future work 


5.5. Comparison of the computation times 

Table 1 shows a comparison of the computation times needed 
to simulate the scenario of Fig. 5 when using the knowledge-based 
model and the black-box model. The knowledge-based simulator 


In this contribution, a general framework for the identification 
of the dynamic of the key variables in a thermodynamic engine is 
proposed. This enables a fast simulator of the engine to be derived 
which can be incorporated in the simulator of the complete hybrid 
power generation station for simulation and/or control design and 
validation purposes. The data used for the identification may come 


Q (rad/sec) J (mm) Tu (K) 



Fig. 9. Profil of the exogenous signal that is used in the validation step. 


true and identified work true and identified heat exchanged 




Fig. 10. validation data: time profiles of the true and the predicted work and heat exchange. 
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either from purely experimental measurement or from knowledge- 
based model that one would like to match in order to accelerate the 
simulation and the control design step. 

Although the validation of the proposed solution has been 
shown in the case where the variables to be identified are the 
work and the efficiency, any other interesting variables charac¬ 
terizing the cycle can be identified such as the excursion of the 
pressure/temperature or the maximum values of these variables. 
The latter can be crucial in order to meet the constraints-related 
requirement. 


Appendix A. 

A. 1. Definition of the basis functions B^(rj) 


The basis function involved in the definition of the map r are 
given by [3]: 


:= {l}u<cr 


(O 


n m -l 


,(0 


n m 


(A.1) 


where the number of functions in the set if equal to n& = 2 n m while 
the functions and B^ are given by: 


:=(1 +«i)j 


V 


+ a 


BfW ■= 


rj 


1 +£*,-(1 - r ]) 


(A.2) 


The coefficients are given by a x ■: = - 1 for some constant 

0. 
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